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^ Abstract 

We report here on the recent application of a now classical general reduction 

5—1 technique, the Reduced-Basis (RB) approach initiated in [PRV+02 , to the spe- 

ed 



cific context of differential equations with random coefficients. After an elemen- 
tary presentation of the approach, we review two contributions of the authors: 
[BBM+09 , which presents the application of the RB approach for the discretiza- 



tion of a simple second order elliptic equation supplied with a random boundary 
condition, and [BL09J, which uses a RB type approach to reduce the variance 
in the Monte-Carlo simulation of a stochastic differential equation. We con- 
clude the review with some general comments and also discuss possible tracks 
for further research in the direction. 
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1 Introduction 



In this work we describe reduced basis (RB) approximation and a posteriori 
error estimation methods for rapid and rehable evaluation of input-output rela- 
tionships in which the output is expressed as a functional of a field variable that 
is the solution of an input-parametrized system. In this paper our emphasis is 
on stochastic phenomena: the parameter is random; the system is a partial dif- 
ferential equation with random coefficients, or a stochastic differential equation, 
namely a differential equation forced by a Brownian process. 

The reduced basis approach is designed to serve two important, ubiquitous, 
and challenging engineering contexts: real-time, such as estimation and control; 
and many-query, such as design, multi-scale simulation, and — our emphasis 
here — statistical analysis. The parametric real-time and many-query contexts 
represent not only computational challenges, but also computational opportu- 
nities: we may restrict our attention to a manifold of solutions, which can be 
rather accurately represented by a low- dimensional vector space; we can accept 
greatly increased pre-processing or "Offiine" cost in exchange for greatly de- 
creased "Online" cost for each new input-output evaluation. (All of these terms, 
such as "Online," will be more precisely defined in the Section 2.1 which consti- 
tutes a pedagogical introduction to the reduced basis approach.) Most variants 
of the reduced basis approach exploit these opportunities in some important 
fashion. 

Early work on the reduced basis method focused on deterministic alge- 
braic and differential systems arising in specific domains |FM71[ IASB78[ INoo81[ 
INoo82[ INPSOi INP83a[ INP83bl |Nag79| ; the techniques were subsequently ex- 
tended to more general finite-dimensional systems as well as certain classes of 
partial differential equations (and ordinary differential equations) |BR95[ IFR831 
ILinQll INoor84i INPA84[ IPorSTi IRhe81[ IRheQSl IPor85) : the next decades saw 
further expansion into different applications and classes of equations, such as 
fluid dynamics and the incompressible Navier-Stokes equations |Pet89[ IGun89] . 
There is ample evidence of potential and realized success. 

Recent research in reduced basis methods for deterministic parametrized 
partial differential equations both borrows from earlier efforts and also empha- 
sizes new components: sampling techniques for construction of optimal reduced 
basis approximation spaces in particular in higher dimensional parameter do- 
mams [SenOSI lBBM+n9[ |Nguyen07| ; rigorous a posteriori error estimation in 
appropriate norms and for particular scalar outputs of interest |KP09|, IHOOSa] ; 
and fastidious separation between the offline stage and online stage of the com- 
putations to achieve very rapid response |N VP05| . These reduced basis methods 
can now be applied to larger, more global parameter domains, with much greater 
certainty and error control. 

In this paper we emphasize the application of certified reduced basis meth- 
ods to stochastic problems. Two illustrative approaches are explored. In the 
first approach [BBM+09j we consider application of the reduced basis approach 
to partial differential equations with random coefficients: we associate realiza- 
tions of the random solution field to deterministic solutions of a parametrized 
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deterministic partial differential equation; we apply the classical reduced basis 
approach to the parametrized deterministic partial differential equation. Sta- 
tistical information may finally be obtained, for example through Monte Carlo 
approximations. New issues arise related to the simultaneous approximation of 
both the input random field and the solution random field. 

In the second approach |BL09| we directly consider a statistical embodiment 
of the reduced basis notions. Here reduced basis ideas originally conceived in the 
deterministic differential context are re-interpreted in the statistical context: the 
deterministic differential equation is replaced by a parametrized random process; 
snapshots on the parametric manifold are replaced by correlated ensembles on 
the parametric manifold; error minimization (in the Galerkin sense) is replaced 
by variance reduction; offline and online stages are effected through fine and 
coarse ensembles. This technique is here applied to parametrized stochastic 
differential equations. 

We begin, in Section [2] with an initiation to the RB approach, considering a 
simple, prototypical elliptic problem, with deterministic coefficients. Section [3] 
then presents the application of the approach to a boundary value problem sup- 
plied with a random boundary condition. The section summarizes the results 
some of us obtained in [BBM+09| . With Section [i) we address a problem dif- 
ferent in nature, although also involving randomness. The issue considered is 
the variance reduction of a Monte-Carlo method for solving a stochastic differ- 
ential equation. The RB approach has been successfully employed in |BL09) 
to efficiently generate companion variables that are used as control variate and 
eventually reduce the variance of the original quantities. The section outlines 
the approach and shows its success on representative results obtained. We con- 
clude the article presenting in Section [5] some potential, alternate applications 
of the approach in the random context. 

2 An Initiation to Reduced-Basis Techniques 

We begin with an overview of Reduced Basis techniques. The level of our expo- 
sition is elementary. Our purpose here is to introduce the main ideas underlying 
the approach, leaving aside all unnecessary technicalities. The reader already 
familiar with this family of approximation approaches may easily skip this sec- 
tion and directly proceed to sections[3]and|4]where the adaptation of the general 
technique to the specific case of partial differential equations with random co- 
efficients and to variance reduction using the RB approach will be addressed. 
We also refer to |RHP[ |Qua09| for pedagogic introductions to the standard RB 
method, though with different perspectives. 

2.1 Outline of the Reduced Basis approach 

Assume that we need to evaluate, for many values of the parameter /Lt, some out- 
put quantity s(/x) = F{u(fi)) function of the solution u{fi) to a partial differential 
equation parametrized by this parameter /i. If the computation of and s(/i) 
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for each single value of the parameter n already invokes elaborate algorithms, 
and this is indeed the case in the context of partial differential equations, then 
the numerical simulation of and s(^) for many /j, may become a computa- 
tionally overwhelming task. Reducing the cost of parametrized computations is 
thus a challenge to the numerical simulation. This is the purpose of Reduced 
Basis technique (abbreviated as RB throughout this article) to reduce this cost. 

Let us formalize our discussion in the simple case of a partial differential 
equation which is an elliptic second order equation of the form (see ([6| below): 



on a domain V with homogeneous Dirichlet boundary conditions. The mathe- 
matical setting is classical. We assume that the solution of interest u{fi) € X is 
an element of a Hilbert space X with inner product (•, •)x and norm || • 
The output s(/i) = F{u{fi)) e M is a scalar quantity where F : X ^ R 
is a smooth (typically linear) function and ^ is a P-dimensional parameter 
varying in a fixed given range A C M^. An example of such output s is 

s(/i) = F{u{fi)) :— / f u{ii) (see (8) below). The function u{^) is mathe- 

matically defined as the solution to the general variational formulation: 

Find ii(/i) G X solution to a(u(/i), v; ji) — l{v) , Vu G A" , (1) 

where a(-, •; /i) is a symmetric bilinear form, continuous and coercive on X and 
where is a linear form, continuous on X. For all e A, a(-, •; /i) thus defines 
an inner product in X. The existence and uniqueness of ■ti(/i), for each /x, is 
then obtained by standard arguments. 

We henceforth denote by || • ||^ the norm || • jj^ = \/a[-, •; ^) equivalent to || • ||x 
(under appropriate assumptions on A, see below), which is usually termed the 

energy norm. In the sequel, we denote by uj\f{^) e Xj^ an accurate Galerkin 
approximation for u(/i) in a linear subspace Xj^ C AT of dimension A/" S> 1 
and by s_\f{fi) = F{u^f{fi)) the corresponding approximation for the output 
s{fi). For that particular choice of ATjv, we assume that the approximation 
error |s(/i) — s_\f{p) \ is uniformly sufficient small for all /i G A. That is, s_\f{p) is 
considered as a good approximation of the output s(/i) in practical applications. 
The difffculty is, we put ourselves in the situation where computing sj\f{ii) for 
all the values fi needed is too expensive, given the high dimensionality M of the 
space X_\f and the number of parameters fi for which equation ([ij need to be 
solved. 

The RB approach typically consists of two steps. The purpose of the first 
step is to construct a linear subspace 

X^,jv = Span(uAA(M^),n = l,...,Ar) , (2) 

subset of X_\f, of dimension N ^ J\f, using a few approximate solutions to ([I]) 
for particular values of the parameter /i. The point is of course to carefully 
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select these values {^J'n)l<n<N S A''^ of the parameter, and we will discuss this 
below (see Q and For intuitively clear reasons, the particular solutions 

u^{fjL^) are called snapshots. This first step is called the offline step, and is 
typically an expensive computation, performed once for all. In a second step, 
called the online step, an approximation uj^^^if^) S Xn',n of the solution to ([T]) 
is computed as a linear combination of the ujs/ [fi^ ) . The problem solved states: 

Find uj^,n{^i) e Xjs/^n solution to a{u_sf^pf{^), v; fi) — l{v) , Vv G Xj^^n ■ (3) 

This problem is much less computationally demanding than solving for the fine 
solution u_\f(fi), and will be performed for many values of the parameter /x. 
We denote by sj^^n{^) — F{uj\f^NilJ')) the corresponding approximation of the 
output s(/i). An a posteriori estimator A^(/i) for the output approximation 
error \sjij-{iJ,) — s_\f,N{^j,)\ is needed in order to appropriately calibrate N and 
select the {fJ.n)i<n<N ■ This a posteriori estimator may also be used in the 
online step to check the accuracy of the output. We shall make this precise 
below. For the time being, we only emphasize that the a posteriori analysis we 
develop aims at assessing the quality of the approximation of the output s (and 
not the overall quality of the approximation of the solution), see |AO00) and 
references therein. The method is typically called a goal oriented approximation 
method. 

The formal argument that gives hope to construct an accurate approxi- 
mation of the solution u{ii) to ([T]) using this process is that the manifold 
A^jv" — {^^aa(/^)i A* € A} is expected to be well approximated by a linear space 
of dimension much smaller than A/", the dimension of the ambient space Xj\j-. 
An expansion on a few snapshots N has therefore a chance to succeed in ac- 
curately capturing the solution u(fj,) for all parameter values fi. The reduced 
basis method is fundamentally a discretization method to appoximate the state 
space Ai^f, with a view to computing an accurate approximation of the output. 
Of course, this requirement strongly depends on the choice of the parametriza- 
tion which is a matter of modelling. 

The RB method yields good approximations SAA.Ar(/i) of s_\f{fi) under appro- 
priate assumptions on the dependency of the solution u(/i) on the input param- 
eter fi. As a consequence, optimal choices for the approximation space Xjvjv 
should account for the dependency of the problem with respect to /i. More pre- 
cisely, the method should select parameter values {Hn)i<n<N G A^ with a view 
to controlling the norm of the output approximation error |s7v(/i) — sj\f^N{fi)\ 
as a function of /i. For most applications, the appropriate norm to consider for 
the error as a function of /i is the L°° norm and this is the choice indeed made 
by the RB approach, in contrast to many other, alternative approaches. The 
desirable choice of (/i^)i<n<7V is thus defined by: 

(M«)i<ri<we arginf I sup |saa(m) ~ ^^.^(m)! I (4) 

Note that, although not explicitly stated, the rightmost term saA.A'(m) in Q 
parametrically depends on {tJ-n)i<n<N because the solution to ([T]) for fi is de- 
velopped as a linear combination of the corresponding snapshots u^{fin). 
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It is unfortunately very difficult to compute Q in practice. With the publi- 
cation [PRV+Oi] . the RB approach suggests an alternative, practically feasible 
procedure. Instead of the parameters {fJ-n)i<n<N defined by Q, the idea is to 
select approximate minimizers of 

{^J'n)l<n<N e arginf ( sup A%{fi) ) . (5) 

Note that there are two differences between Q and ([S]). First, the set A has 
been discretized into a very large trial sample of parameters Atrial C A. Second, 
and more importantly, the quantity A^(/x) minimized in ^ is an estimator of 
IsN'ifJ') ~ ^N{fi)\. A fundamental additional ingredient is that the approximate 
mimmizers of ^ are selected using a specific procedure, called greedy because 
the parameter values /i^, n = 1,...,N, are selected incrementally. Such an 
incremental procedure is in particular interesting when N is not known in ad- 
vance, since the computation of approximate pi^ (1 < « < N) does not depend 
on N and may be performed until the infimum in ^ is judged sufficiently low. 

Of course, the computation of approximations to ([5| with such a greedy al- 
gorithm can still be expensive, because a very large trial sample of parameters 
Atrial C A might have to be explored. The RB method is thus only considered 
efficient when the original problem, problem ([ij here, has to be computed for 
such a large number of input parameter values /i, that the overall procedure 
(computationally expensive offline step and then, efficient online step) is prac- 
tically more amenable than following the original, direct approach. One often 
speaks of a many-query computational context when it is the case. Notice that 
the RB method is not to be seen as a competitor to the usual discretization 
methods; it rather builds upon already efficient discretization methods using 
appropriate choices of X_\f in order to speed up computations that have to be 
performed repeatedly. 

The approach can therefore be reformulated as the following two-step pro- 
cedure 

• in the offline stage (which, we recall, may possibly be computationally 
expensive), one "learns" from a very large trial sample of parameters 
Atrial C A how to choosc a small number N of parameter values; this 
is performed using a greedy algorithm that incrementally selects the 

n = 1, . . . , A^; the selection is based on the estimator A^(/x); accurate ap- 
proximations u_\f[fin) for solutions w(/i„) to ([T]) are correspondingly com- 
puted at those few parameter values; 

• in the online stage, computationally inexpensive approximations UA^,Ar(/x) 
of solutions w(/Lt) to ([T]) are computed for many values /i G A of the param- 
eter, using the Galerkin projection Q; the latter values need not be in 
the sample Atrial, and yield approximations s^^n{^) for the output s(/i); 
the estimator A|^(/i), already useful in the offline step, is again employed 
to check the quality of the online approximation (this check is called cer- 
tification). 
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Notice that the computation of the error estimator (^) needs to be inexpen- 
sive, in order to be efficiently used on the very large trial sample in the offline 
stage, and for each new parameter values in the online stage. 

One might ask why we proceed with the reduced basis approach rather than 
simply interpolate s(/i), given the few values {s_\f{^i), . . . , s_\f{^N)}. There are 
several important reasons: first, we have rigorous error estimators based on 
the residual that are simply not possible based on direct interpolation; second, 
these residuals and error bounds drive the greedy procedure; third, the state- 
space approach provides Galerkin projection as an "optimal" interpolant for the 
particular problem of interest; and fourth, in higher parameter dimensions (say 
of the order of 10 parameters), in fact the a priori construction of scattered-data 
interpolation points and procedures is very difficult, and the combination of the 
greedy and Galerkin is much more effective. 

We are now in position to give some more details on both the offline and 
online steps of the RB approach in a very simple case: an elliptic problem, with 
an affine dependency on the parameter. Our next section will make specific what 
the greedy algorithm, the estimator A^(/i), along with other objects abstractly 
manipulated above, are. 



2.2 Some more details on a simple case 

As mentioned above, we consider for simplicity the Dirichlet problem 

-div(4(Ai)V«(M)) =/inI?, 

V / (6) 

u{fi) = on dV , 

where 2? is a two-, or three-dimensional domain and the matrix A{p) is param- 
eterized by a single scalar parameter /x € A = [^J-mim fJ-max] C M*^. We assume 
that the matrix A is symmetric and depends on ^ in an affine way: 

Aif^) ^ Ag + fiA^, Vfie A. (7) 

This assumption ([T]) is a crucial ingredient, responsible, as we shall explain 
below, for a considerable speed-up and thus for the success of the RB approach 
here. More generally, either we must identify by inspection or construction an 
"affine" decomposition of the form ([t]) , or we must develop an appropriate affine 
approximation; both issues are discussed further below. 

We assume we are interested in efficiently computing, for many values of 
fi € A, the output: 

s{fi) = F{u{fi)) :- / /u(^). (8) 
Jv 

This is of course only a specific situation. The output function can be much 
more general, like a linear form / gu{^) with some g f- Many other cases 
are possible, but they all come at a cost, both in terms of analysis and in terms 
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of workload. The case (|8]), where the output coincides with the hnear form 
present in the right-hand side of the variational formulation of ^ (and where 
the bilinear form a involved in the variational formulation is symmetric), is 
called compliant. Having ([S]) as an output function in particular simplifies the 
a posteriori error analysis of the problem (namely the construction of A^(/i)). 

We equip the problem, somewhat vaguely formulated in ([6])-([7])-(|8]) above, 
with the appropriate mathematical setting that allow for all our necessary ma- 
nipulations below to make sense. For consistency, we now briefly summarize this 
setting. The domain V is an open bounded connected domain with Lipschitz 
boundary dV, the right-hand side / € L'^i'D) belongs to the Lebesgue space 
of square integrable functions, A{fj,) is a symmetric matrix, which is positive- 
definite almost everywhere in D. Each entry of A{fi) is assumed in L°°{'D). We 
assume Aq is symmetric positive-definite, and Ai is symmetric positive. The 

ambient Hilbert space X is chosen equal to the usual Sobolev space Hq (2?) . The 
function u(/i) is defined as the solution to the variational formulation ([T]) with 

a{w, v\ jjL) = I A(/x)Vw • Vw, l{v) = f V, for all v, w, in X and all /i G A. 
Jv Jv 
As for the direct discretization of the problem, we also put ourselves in a 

classical situation. If V is polygonal for instance, there exist many discretiza- 
tion methods that allow to compute Galerkin approximations uj\f(fi) of u(/i) 
in finite dimensional linear subspaces Xj^/ of X for any fixed parameter value 
/z G A. The Finite-Element method |SF73| ICia78| is of course a good exam- 
ple. Then, for each parameter value fj, G A, the numerical computation of 
uj\f{^) — Yl^=iUn{lj)4>n on the Galerkin basis {4>n)i<n<N' of Xj^ is achieved 
by solving a large linear system 

Find U{^i) e R-^ solution to E{fi)U{fi) =6, 

for the vector [/(/i) = iUn{p))i<n<N' ^ ^ where b = il{(f>n))i<n<Ar a vector 
in M.-^ and B{p) = + /i -Bi is a A/" x A/" real invertible matrix. Note that 

the assumption of afhnc parametrization Q makes possible, for each parameter 
value /i, the computation of the entries of the matrix B_{^) in 0{M) operations 

(due to sparsity), using the precomputed integrals {Bq).. = J^AqWcj)i ■ '^(j>j, 

i,j = 1, . . . ,J\f for q = 0,1. The evaluation of U{^) for many J ^ 1 parameter 
values yU using iterative solvers costs J x 0{J\f'^) operations with fc < 3 jGvL96] . 
where k depends on the sparsity and the conditioning number of the involved 
matrices. 

As mentioned above in our general, formal presentation, the goal of the RB 
approach is to build a smaller finite dimensional approximation space Xf^f^^ C 
Xj^ sufficiently good for all G A, with ^ A^, so that the computational cost 
is approximately reduced to iV x 0{N^) + J x 0{N^), where N x 0{N^) is the 
cost of offline computations and J x 0{N^), the cost of online computations, is 
independent of M , using the Galerkin approximation ([S]) in X_\f,N. 
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We now successively describe in the following three paragraphs the construc- 
tion of the a posteriori estimator, that of the greedy algorithm employed in the 
offline step, and the combination of all ingredients in the online step. 



2.2.1 A posteriori estimator. 

For the coercive elliptic problem (|6]), the a posteriori error estimator A^(/^) for 
the output RB approximation error \sf^{ii) — sj\f_N(fi) \ is simple to devise, based 
on a global a posteriori error estimator Ajv(/^) for ||M7v(/i) — u^^,N{^J■)\\^J, using a 
classical technique with residuals. 

We refer to VRPOl IVPHPOS^ [HPOTl |Dep09[ IPEOTal INHHP091 INHPOSj 
for the construction of similar a posteriori error estimators in various applied 
settings of the RB method. 

We first define the residual bilinear form 

g{w, v; /i) = a{w, v; — l{v) , Vw, u € X , V/i € A , 

and the operator G{^) : Xj^ Xjs/ such that 

g{w, v; fi) = (G(^) w, v)^ , Vw, w G , G A . 

We next assume we are given, for all /i G A, a lower bound aLsilJ) for the 
coercivity constant of a(-, •; /i) on A"jv, that is, 

< aLB(M) < «c(m) = inf ^tS^ , G a . (9) 

The lower bound ai^sif-) can be given by an a priori analysis before discretiza- 
tion {otLBip-) would then be the coercivity constant of a(-, •; /i) on X), or numer- 
ically evaluated based on an approximation procedure, which might be difficult 
in some cases, see |HRSP07I IRHP] . 

Then the a posteriori estimator we use is defined in the following. 

Proposition 2.1. For any linear subspace X_\f^]y of X_\f, there exists a com- 
putable error bound A^(/l() such that: 

Mt^) - smM < A^(A^) := , Vm G A . (10) 

For consistency, we now briefly outline the proof of this proposition. We 
simply observe the sequence of equalities 

|s7v(m) - stV.jvIm)! = \F{uAf{fJ-)) - ^("aa,w(m))I 

^\\u^f{^J-)-u^f,N{^J■)\\l■ (11) 
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using the linearity of F — Z, the variational problem and its discretized ap- 
proximation in X^, the symmetry of a(-,-;/i) and the fact that a{uj^{^) — 
'Uj\f,N{fJ-),v) = 0, for all v G Xj^/^pf. On the other hand, inserting v — uj\f{fj,) — 
uj^,N{y) in the general equality a{u^{fi) - u^^N{fi),v; ^) = -g{uj^ ^m{p)-,v] n) 
(for all V S Xj^), and using the bound ^/aijj{ji)\\v\\x < (for aU v e Xj^), 
we note that 

II ( \ ( \\\ ^ A I \ I|G'(m) UAA,Ar(M)IU no^ 
||uaa(m) - UM,N{m^i < AAr(^) := , ' , ■ (12) 



We conclude the proof of combining (11) with (12). 



We may similarly prove (but we will omit the argument here for brevity) the 
inverse inequality: 

a^(m) < f |.saa(m) - s^M^^)\ , (13) 

using the continuity constant 

IKN = sup sup II |, II |, (14) 

ioeXj^\{Q}vex,v\{o} \\w\\x\\v\\x 

of the bilinear form a(-, •; fi) on Xj\f for all /i € A, which is bounded above by 



the continuity constant on X. The inequality (13) ensures sharpness of the a 
posteriori estimator ( |10[ ) , depending of course on the quality of the lower-bound 

2.2.2 Offline stage and greedy algorithm. 

The greedy algorithm employed to select the snapshots u^{fin) typically reads: 



1: choose /xi e A randomly 

2: compute uj\f{iJ,i) to define Xj\f,i = Span (uAr(/ii)) 
3: for n = 2 to do 

3: choose fin S argmax{ A^_i (^) , ^ e Atrial} 

3: compute uj\f{fin) to define Xj^/^n = Span {uj\f{fim) , m — 1, . . . ,n) 
4: end for 

In the initialization step, we may equally use /J,i £ argmax{|s(/i)| , /i G Agniaiitriai}) 
where Agmaiitriai C A is a very small trial sample in A, much smaller than A itself. 
Likewise, the algorithm can in practice be terminated when the output approx- 
imation error is judged sufficiently small (say, |A^(/z)| < e for all e Atrial), 
and not when the iteration number reaches a maximum n = N . 

The choice of the trial sample Atrial (and similarly, the smaller sample 
Asmaiitriai) IS a dclicate practical issue. It is often simply taken as a random 
sample in A. Of course, this first guess may be insufficient to reach the required 
accuracy level e in A^(/x), for all G A, in the online stage. But fortunately. 
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if the computation of A^(/i) for any S A is sufficiently inexpensive, one can 
check this accuracy onHne for each query in jj,. Should A^(/i) > e occur for 
some online value of the parameter fi, one can still explicitly compute uj\f{^) for 
that exact same n and enrich the space Xj^jy correspondingly. This bootstrap 
approach of course allows to reach the required accuracy level e at that /it. It pro- 
vides significant computational reductions in the online stage provided that the 
RB approximation space Xj^f^iq does not need to be enriched too often online. 
We will explain the methodology for fast computations of A|^,(/^) below. 

The offline selection procedure needs to be consistent with the online proce- 
dure, and thus the above greedy algorithm uses the same estimator A^(/i) for 
all /i S Atrial as the online procedure. Since the computation of A^(/Lt) is, by 
construction and on purpose, fast for all /i G A, the exploration of a very large 
training sample Atrial (which is a subset of A) is possible offline. 

No systematic procedure seems to be available, which allows to build good 
initial guesses Atrial e^; nihilo. Even for a specific problem, we are not aware 
either of any a priori results that quantify how good an initial guess Atrial is . 
The only option is, as is indeed performed by the RB approach, to a posteriori 
check, and possibly improve, the quality of the initial guess Atrial (however, 
the quality of the initial guess Atrial can be slightly improved offiine by using 
adaptive training samples in the greedy algorithm ^HOOSbj ). 

The estimators A^(/i) are employed in the greedy algorithm to filter can- 
didate values for A. Numerous numerical evidences support the success of this 
pragmatic approach |VRPn2L IVPHP081 iHPOTl IPH07bl IPROTal |Dep09[ INHP08L 
INRHPOQj . 

Last, notice that the cost of offline computations scales as Woffiinc = 0(|Atriai|) x 
(j2n=i ^online ('T')^ -|- X 0{N^) where Woniine('T') is the marginal cost of one 
online- type computation for UJ^^n{,^^) and A^(/i) at a selected parameter value 
A* G Atrial (where 1 < n < — 1), and 0(| Atrial |) includes a max-search in Atrial- 
(Recall that fc < 3 depends on the solver used for large sparse linear systems.) 

2.2.3 Online stage: fast computations including a posteriori estima- 
tors. 

We now explain how to efficiently compute uj^^nip), s^f^n{^J■) and Afj(/i) once 
the RB approximation space X_\f_n has been constructed. This task has to 
be completed twice in the RB approach. First, this is used in the many of- 
fline computations when e Atrial explores the trial sample in order to find 
yU„ S argmax{A^_]^(yLt), /X €: Atrial} at each iteration n of the greedy algorithm. 
Second, this is used for the many online computations (when n = N). We 
present the procedure in the latter case. 

By construction, the family (wA/'(Mn))i<n<Af generated by the greedy algo- 
rithm described in Section r2. 2. 21 is a basis of 

Xatm = Span {uj^(p.n) , n = l,...,N) . 
For any /i € A, we would then like to compute the RB approximation Uf^^^di) = 
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"^^=1 UN,n{fJ')'U'M{l^n), which can be achieved by solving a small N x N (full) 
linear system 

Q{li)UN{ij) = c, 

for the vector Uniy) = (f^w,"(M))i<„<7v ^ ^ ^'^^^ ^ = (K"AA(A*n)))i<„<Ar a 
vector in and CZ(a') = Co+^Ci is a A^x real invertible matrix. In practice, 

the matrix C{^) is close to a singular matrix, and it is essential to compute the 

RB approximation as u^^n{^) = YlH^i UN.n{y)Cn using a basis (Cn)i<n<Ar of 
Xf^^N that is orthonormal for the inner-product (•, ■)x- The determination of 
appropriate {Cn)i<n<N is easily performed, since N is small, using Simple or 
Modified Gram-Schmidt procedures. The problem to solve states: 

Find UNif^) e solution to Q{fi)UN{i^) = 5, (15) 
where C/7v(a') = {UN.7iif^)) e K^, c (KCri))i<„<Ar is a vector in 

\ / l<n<N — — 

and C(/x) = Cq_ + fJ.Cn is a N x N real invertible matrix. So, for each parameter 

value /I e A, the entries of the latter matrix £Z(/i) can be computed in 0{N'^) 

operations using the precomputed integrals {Gq)^. — j^^Aq^Qi ■ VCj, i,j = 

\, . . . ,N for q — 0, 1. (Note that the assumption of affine parametrization is 
essential here). And the evaluation of U]\[{^) for many J ^ 1 parameter values 
H G A finally costs J x 0{N^) operations using exact solvers for symmetric 
problems like Cholesky |GvL96j . 

For each /i e A, the output sj\f^]s[{^) = F(uj\f^N{ijL)) can also be computed 
very fast in 0{N) operations upon noting that F is linear and all the values 
Fi^M ,N{fJn)), n = 1,. . ■ ,N can be precomputed offline. The corresponding a 



posteriori estimator A^(/x) given by (10) has now to be computed, hopefully 
equally fast. Because of the affine dependence of ^(/i) on fi, a similar affine 

dependence G{fj,) = Ga + l-J. Gi holds for the operator G (for all fj, € A), where 

(Go w,v)x — / AqVw ■ Vu — / fv, and (Gi w,v)x = / ^iVw • Vw for all 

Jv Jv Jv 

!>, w, in Xj\f. So one can evaluate very fast the norm 

2 



||G(^) Uf^^N{fl)\\x = \\Go WA/-,Ar(/i)||x + 2^(Go U^f,N{fj), Gi u^f^xip)) 



X 



+ M'l|Giu^,jv(/i)llx (16) 

for /I e A, once, with obvious notation, the scalar products (GiUf^^NifJ-p), GjUj^^^^f. 
have been precomputed offline and stored. Assuming that the lower-bound 



Q^LB (m) used in ( 10 1 is known, the computation of the a posteriori estimator 
A|^(/i) itself is thus also very fast. Notice that the affine parametrization ([7| 
plays a crucial role in the above decomposition of the computation. 

Finally, the marginal cost of one online-type computation on X_\f_n for one 
parameter value ^ is WonVmc{n) — 0{v?) (where n — 1, . . . , A^). So, assuming 
that no basis enrichment is necessary during the online stage using the RB 
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approximation space Xj^^n (that is, Aff{p) < e for all the parameter values /x 
queried online) , the total online cost for many J ^ 1 parameter values fi scales 
as Woniinc ~ J X 0{N^). And, the total cost of computations with the RB 
approach is then Woffiinc + T^o„ii„c = N x 0{N^) + [J + 0(|AtHai|)) x 0{N^), 
which has to be compared to J x 0(N^) operations for a direct approach (with 
fc < 3 depending on the solver used for large sparse linear systems) . In the limit 
of infinitely many online evaluations J 3> 1, the computational saving of the 
RB approach is tremendous. 



2.3 Some elements of analysis, and some extensions, of 
the RB method 

2.3.1 Some elements of theory. 

The RB approach has undoubtedly proved successful in a large variety of ap- 
plications IVRP021 IVPRP031 iHPOTl IPROTbl IPROTal |Dep09[ INRP081 INRHPOQj . 

The theoretical understanding of the approach is however still limited, and is 
far from covering all practical situations of interest. Of course, little theory is to 
be expected in the usual a priori way. As already explained, the RB approach 
is deliberately designed to a posteriori adapt to practical settings. The only 
available a priori analysis is related to two issues: the expected "theoretical" 
quality of the RB approximation, and the efficiency of the greedy algorithm. 
We now briefly summarize what is known to date on both issues. 

The RB approach is in fact expected to perform ideally^ in the following 
sense. In the context of our simple problem ([6]) , it is possible, adapting the clas- 
sical Lagrange interpolation theory to the context of parameterized boundary 
value problems and assuming that the matrix Ai is non-negative, to obtain an 

upper bound of Q. The following theoretical a priori analysis result follows. 
It states the exponential accuracy of the RB approximation in terms on the 
dimension N of the reduced basis. 

Proposition 2.2. For all parameter ranges A :~ [Mmin, Mmax] C M^ , there 
exists an integer Nq = O (in ^^7^^^ as — > +00, and a constant c > 
independent of A such that, for all N > Nq > 2, there exist N parameter values 
Mmin —■ < . . . < < X^^i < . . . < := /imaxj n — 2, . . . , N — 2, 

sastisfying (recall \\ ■ ||o = || • ||^ with /i = is an Hilbertian norm on X): 

sup (inf {||uA-(/i) - w\\o , w e Span {u^r{\'^),n = 1,...,N)}) 

<e-«^(^-i)sup||u^(/i)||o. (17) 



We refer to |Boy09[ Chapter 4] and jPR07b| for the proof of Proposition [23 

The approximation space Span (^u_\f{X^),n = 1, . . . ,N) used for the state- 
ment and the proof of Proposition |2.2| is different from the RB approximation 
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space X^f^N built in practice by the RB greedy algorithm. Numerical exper- 
iments even suggest that it is not an equally good choice (see |PR07bp . So 
it is desirable to better understand the actual outcome of the RB greedy al- 
gorithm used offline. The concept of greedy algorithm appears in many nu- 
merical approaches for problems of approximation. It typically consists in a 
recursive procedure approximating an optimal solution to a complex problem, 
using a sequence of sub-optimal solutions incrementally improved. Otherwise 
stated, each iteration takes the solution of the previous iteration as an initial 
guess and improves it. In the theory of approximation of functions in par- 
ticular |Dev931 IV.N08| , greedy algorithms are used to incrementally compute 
the combinations of functions from a given dictionnary which best approximate 
some given function. The RB greedy algorithm has a somewhat different view- 
point: it incrementally computes for integers N some basis functions uj\f{iJ,n), 
n = 1, ... ,7V, spanning a linear space Xj^/^n that best approximates a family 
of functions uj^{fi), £ A. The RB greedy algorithm however has a flavour 
similar to other greedy algorithms that typically build best-approximants in 
general classes of functions. It is therefore possible to better understand the 
RB greedy algorithm using classical ingredients of approximation theory. The 
notion of Kolmogorov width is an instance of such a classical ingredient. We 
refer to |Boy09[ Chapter 3] and |BMPPT10| for more details and some elements 
of analysis of the RB greedy algorithm. 



2.3.2 Extensions of the approach to cases more general than (|6|. 

The RB approach of course does not only apply to simple situations like 
Many more general situations may be addressed, the major limitation to the 
genericity of the approach being the need for constructing fast computable a 
posteriori error estimators. 

Instances of problems where the RB approach has been successfully tested 
are the following: affine formulations, non-coercive linear elliptic problems, non- 
compliant linear elliptic problems, problems with non-affine parameters, non- 
linear elliptic problems, semi-discretized (nonlinear) parabolic problems. The 
purpose of this section is to briefly review these extensions of our above simple 
setting. In the next section, we will then introduce a problem with random co- 
efficients. For simplicity, we take it almost as simple as the above problem ([6]), 
see (23)-(24| below. We anticipate that, if they involve a random component, 
most of the extensions outlined in the present section could also, in principle, 
be treated using the RB approach. 



Affine formulations. Beyond the simple case presented above in Section [Z2l 
which involves an elliptic operator in divergence form afflnely depending on the 
parameter, the RB approach can be extended to general elliptic problems with 
variational formulation of the form 

Find m(^) e X solution to g{u{^),v; /i) = , Wv e X , (18) 
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where the form g(-, ■; fi) on X x X admits an affine parametrization, that is, 
writes 

Q 

g{w,v;n) =^^eq{fi)gq{w,v) , yw,v e X , e A , (19) 

9=1 

with parameter-independent forms (ffgl'i ■))i<g<Q (where some of the gq may 
only depend on v) and coefficients (0(;(M))i<q<Q- We emphasize that the whole 
RB algorithm presented in the simple case above directly translates in this 
situation. In particular, the matrices used in the online evaluation procedure 
can be constructed offline. 

Non-coercive symmetric linear elliptic problems. The RB approach can 
be extended to the case where the symmetric continuous bilinear form a(-, ■; fj,) 
is not coercive but only inf-sup stable. An example is the Helmholtz prob- 
lem treated in [SVH+06) . Our discussion of the elliptic problem above can be 
adapted in a straightforward way, the only change in offline and online compu- 
tations being that the inf-sup stability constant on X^: 

< hsip) < /3(m) inf sup ."^"'''j'^^ , V/i G A , (20) 

is substituted for aLBifJ-)- In practice, the evaluation of I3lb{i^) is typically 
more involved than the evaluation of the coercivity constant aLsii^)- We refer 
to HKY+OQ] for an appropriate technique. 



Non-compliant linear elliptic problems. In the particular choices of 
F — I for the output and of symmetric matrices 4(/i) for the definition of 

the bilinear form a(-, •;/j,) correspond to a particular class of problems called, 
we recall, compliant. Non-compliant linear elliptic problems can be treated as 
well, but this is somewhat more technical. These are the cases where, for some 



/Li G A at least, either u(/i) is solution to a weak form (18) with g{v,w;^) = 
a(v,w;iJ,) — l{w), yv,w G X and the bilinear form a{-,-;fi) is not symmetric, 
or the output is s{fi) — F{u{^i)) ^ l{u{fi)) with any linear continuous bmction 
F:X->-R. 

For instance, we explain how to treat the case of a bilinear form a(-, •; fi) that 
is not symmetric, but of course still continuous and inf-sup stable. The analysis 
requires considering the solution to the adjoint problem 

Find ■(/'(m) e X solution to a{v, V'(m); m) = ~F{v) , e X , (21) 

along with the corresponding Galerkin discretization tpAfitJ-) G A_v, the approx- 



imation space X^ j^ft for the solution to (21), and an additional RB approxi- 
mation space X^ C Xj^ of dimension N* <C M . The a posteriori estimator 



obtained is similar to (10), and writes 

'IG(^) u^^n{i^)\\x \\G*{fj,) ipAf,N*{l^)\\x 



pLBili) 

(22) 



15 



where G* is defined from the adjoint problem (21 ) similarly to how G is defined 



from the original problem. Notice that we again used the inf-sup stability con- 



dition ( 20 1 , which indeed holds true after permutation of the arguments v and 
w (since we work in a finite dimensional space Xj^), the value of the inf sup 
constant however being not the same. To build the reduced basis of the primal 
(respectively the dual) problem, in the offline stage, the a posteriori estimator is 
based on \\G{ii) u_\f^]s[{fi)\\x (respectively \\G*{fj,)'ip_\f^N*{l^)\\x)- Apart from the 
above introduction and use of the adjoint problem, the treatment of the non- 
compliant case then basically follows the same lines as that of the compliant 
case. 

Notice that a simple, but less sharp, estimate of the error (namely the left- 

\F{x)\ 



hand side of (22|) can be obtained as sup -j— — \\G{fJ.) uj^ N{^J')\\x■ This 

simple error bound does not involve the solution of any dual problem, and may 
be of interest in particular in the case when multiple outputs are considered. 



However, the primal-dual error bound (22 1 will be much smaller (since it is 
quadratic and not linear in the residual) and in many situations, very easy to 
obtain, since the dual problem is typically simpler to solve than the primal 
problem (it is indeed linear). 

Non-afRne parameters. We have exploited in several places the affine depen- 
dence of A{ii) in ([6| in terms of the coefficient ^. However, there are many cases 

for which the dependency on the parameter is more complicated, as for exam- 
ple, when associated with certain kinds of geometric variations. Extending the 
RB approach to the case of non-affine parametrization is feasible using suitable 
affine approximations. The computation of approximations X]m=i Am (/^)^m(^) 

for functions A{x; /i) (having in mind as an example the prototypical prob- 
lem ([g])), is a general problem of approximation. A possibility, introduced and 
further developed in BNMPOH IGMNP07|. IMNPP09| is to modify the standard 
greedy procedure described above, using interpolation. In short, the approach 

consists in selecting the coefficients ( /3j)f (/i) J of the approximation 

\ /m=l,...,A/ 

1-M[g{-\IJ)] ■■= Em=i/3m (m)5(-;M^) of order M to g{-; {g denoting here a 



general bilinear form, as in ( 18 )-(|19[)) using an interpolation at the so-called 



magic points selected sequentially with 

xi e argmax|5(-;Mf)| , e argmax |5r(-; /x^) - I™_i[g(-; M^)]| , 

for all m = 2, . . . , M. We refer to the contributions cited above for more details. 

Nonlinear elliptic problems For the extension of the RB approach to non- 
linear problems, one major difficulty is again the construction of appropriate 
a posteriori error estimators, which, additionally, need to be computed effi- 
ciently. Several examples of successful extensions are reported on in the lit- 
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erature |VRP02l IVPRP031 iHPOTl |Dep09[ IPROTal INRHP091 INRPOSj . But no 

general theory can of course be developed in the nonlinear context. 



Semi-discretized parabolic problems After time-discretization, parametrized 
parabolic problems can be viewed as a collection of elliptic problems with the 
time variable as an additional parameter. A natural idea is then to build a re- 
duced basis spanned by solutions for given values of the parameter and the time 
variable. Examples of contributions are jGP05[ [CMNP07] . This first approach 
has been improved by techniques combining the RB idea for the parameter with 
a proper orthogonal decomposition (POD) in the time variable, first introduced 
in |HO08a| and further discussed in |KP09| . A route which would be interesting 
to follow could be to try to adapt on-the-fly, as time goes, the reduced basis 
which is the most adapted to the current time. 



3 RB Approach for Boundary Value Problems 
with Stochastic Coefficients 

The first application of the RB approach to a problem with stochastic coefficients 
is introduced in |BBM^09) . The purpose of this section is to overview this 
contribution, in particular showing how the general RB approach needs to be 
adapted to the specificities of the problem. We refer to jBBM+09| for all the 
details omitted below. 



3.1 Setting of the problem 

Let us denote by {ft,J-,¥) a probability space, and by uj £ ft the stochastic 
variable. We consider the stochastic field U{-,uj) that is the almost sure solution 
to 

- div (Aix)VU{x, uj)) ^0, yx eV , (23) 



supplied with a random Robin boundary condition 

n{x) ■ A{x)\7U{x, uj) + B{x, uj) U{x, uj) = g{x) , \fx e dV . (24) 



In (24), the matrix A{x) writes A{x) = a{x)Id where < <j{x) < oo for a.e. 

X € T). Of course, n denotes the outward unit normal at the boundary of 
the smooth domain 2?. The boundary is divided into three non-overlapping 
open subsets: dV = (Fn U Fr U Fb) (see Fig. [l]). The boundary source term 
g is assumed to vanish everywhere except on Fr where it has constant unit 
value: g(x) = lrR,Va; € dV. The scalar random field B{-,uj), parametrizing 
the boundary condition, also vanishes almost everywhere on the boundary dV 
, except on some subset Fb of the boundary dT) with non-zero measure, where 
< femin < B{-,u) < &niax < OO almost surcly and almost everywhere. Note 
that on Fn, (|24| thus reduces to homogeneous Neumann conditions. Physically, 



U{-,uj) models the steady-state temperature field in a heat sink consisting of 
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Figure 1: V has the geometry of a (piece of) heat sink: a spreader I?2 with a 
fin I?i on top. 



an isotropic material of thermal conductivity a, contained in the domain T). 
The sink is subject to zero heat flux on Fn, a constant flux on Fr modeling 
the heat source, and a convective heat transfer on Fb- The Biot number B 
models the effect of the exterior fluid convection on the solid thermal conduction 
problem inside T). In real world engineering applications, the value of B is only 
approximately known. It is therefore legitimate to encode the uncertainties on 
B using a random field B{-,uj), see |LL02j for more details. 



Correspondingly, the solution to (23l-(24), along with any output computed 



from this solution, are also random quantities. Only statitistics on these quanti- 
ties are relevant. We thus consider two statistical outputs for the problem: the 
expected value E(iS') and the variance Var(S') of the random variable 

S{u;)^J^{U{-,u;))^ [ U{- ,uj) (25) 



linearly depending on the trace of the solution C/( • ,a;) on Fr . 

A typical question, example of an Uncertainty Quantification problem, is to 
quantify the sensitivity of the output S (w) . M any existin g contributions already 
addressed the issue: |BTZn51 IDBOOll IMK051 iDNP+oij . 

A possible approach (which we will indeed adopt here) is to evaluate £(5*) 
and Var(S') with the plain Monte-Carlo method using M independent random 
variables {S"^)i<,-n<M with the same distribution law as S. The expectation 
and the variance are respectively approached by the empirical sums 



m— 1 71 — 1 

(26) 

where the normalization factors used (respectively — and ) allow, as is 

traditional in the community of Monte-Carlo methods, to have unbiased esti- 
mators: E(£;m[(S"")]) = E(5) and E (T4/[(5'")]) = Var(S') for aU M. Large 
values of M are typically needed to obtain from \2Q\ accurate approximations of 
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£(5) and Var(5'). Since, for each m = 1, . . . , M, a new realization of the ran- 



dom parameter B is considered and the boundary value problem (23)-(24) has 



to be solved, the task is clearly computationally demanding. It is a many-query 
context, appropriate for the application of the RB approach. 



3.2 Discretization of the problem 

We begin by considering the Karhunen-Loeve (abbreviated as KL) expansion 

K 

(27) 



B{x, uj) - b G(x) + b J2 *fc(x) Ykiuj) 



fe=i 



of the coefficient B{x,uj) (see |Kar46[ iLoeTSl ISTOGj ). In ([27]), JC denotes the 
(possibly infinite) rank of the covariance operator for B(-,uj), which has eigen- 
vectors {^k)i<k<K and eigenvalues {Xk)i<k<K (sorted in decreasing order). The 
random variables (^/c)i<fc<yc are mutually uncorrelated in Lp(rt) with zero mean, 
G is supposed to be normalized fgjjG = 1 and b = J^dP{uj) Jg^B{-,uj) is a 
fixed intensity factor. 



Based on (271, we introduce the deterministic function 

K 



bix,y)^bG{x) + bY,'^kix) 



Vk 



(28) 



fc=i 



defined for almost all x e dV and all y E C M'^, where A** denotes the range 
of the sequence Y = 0^k)i<k<K random variables appearing in (27). Notice 
that B{x,uj) = b{x,y{uj)). 

It is next useful to consider, for any positive integer K < IC, truncated ver- 
sions of the expansions above, and to define, with obvious notation, [/x(',w) as 
the solution to the problem (23)-(24| where B{-,uj) is replaced by the truncated 
KL expansion BK{-,i^) at order K. Similarly, for all y^ £ A^, UK{-',y^) is 
defined as the solution to 



-div 



iv (^Aix)VuKix; y^ 



, VxeV 



n{x) ■ Aix)yuK{x; y^) + 6^(2;, y^)uK(x; y^) = g{x) , Vx G dV, 



(29) 



where bx is the if-truncated sum ( 28 1 



For a given integer K < K., we approximate the random variable S{uj) by 
S'if(w) := {Uk{-,oj)) where Uk{-t^) = uk{-]Y^ {uj)), and the statistical out- 
puts E(S'k) and Var(S'x) by the empirical sums 

^ M M 
m— 1 n— 1 

(30) 

using M independent realizations of the random vector Y^ . In practice, uk{'', Y^) 
is approached using, say, a finite element approximation UK,j\f{ ■ ',Yj^) with 
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Af ^ 1 degrees of freedom. Repeating the task for M realizations of the K- 
dimensional random vector may be overwhelming, and this is where the RB 
approach comes into the picture. We now present the application of the RB 



approach to solve problem ( 29 1 , parametrized by € A" . 

In echo to our presentation of Section [2j note that problem ( 29 ) is affine in 



the input parameter thanks to the KL expansion ( 28 1 of b, which decouples 



the dependence on x and the other variables. To use the RB approach for this 



problem, we consider S' in ( 30 1 as the output of the problem, the parameter 
being (this parameter takes the values F^'^™, m e {1,...,M} being the 
realization number of Y^) and, as will become clear below, the offline stage is 
standard. On the other hand, in the online stage, the a posteriori estimation is 



completed to take into account the truncation error in X in (28). 

Before we turn to this, we emphasize that we have performed above an 
approximation of the coefficient b, since we have truncated its KL expansion. 



The corresponding error should be estimated. In addition, the problem (29 1 



after truncation might be ill-posed, even though the original problem (23)- (24 1 



is well posed. To avoid any corresponding pathological issue, we consider a 



stochastic coefficient b having a KL expansion (28) that is positive for any 



truncation order K (which is a sufficient condition to ensure the well-posedness 



of (23)-(24)), and which converges absolutely a.e. in dV when K ^ IC. For 
this purpose, (i) we require for fc = 1, . . . , /C a uniform bound ||<I'fe||L°°(rB) — "^i 
(ii) we set T\/AfcZfc with independent random variables Zj. uniformly 

distributed in the range (— -s/S, VS), T being a positive coefficient, and (iii) 
we also ask X^fc^iV-^fc < oo. Note that, if /C = oo, condition (iii) imposes a 
sufficiently fast decay of the eigenvalues while k increases. We will see in 
Section [3^ that this fast decay is also important for the practical success of our 
RB approach. Of course, (i)-(ii)-(iii) are arbitrary conditions that we impose 
for simplicity. Alternative settings are possible. 

3.3 Reduced-Basis ingredients 

We know from Section[2]that two essential ingredients in the RB method are an a 
posteriori estimator and a greedy selection procedure. Like in most applications 
of the RB method, both ingredients have to be adapted to the specificities of 
the present context. 



As mentioned above, the statistical outputs (30) require new a posteriori 
estimators. Moreover, the statistical outputs can only be computed after M 
queries Y^ ^ m = 1,...,M, in the parameter y^, so these new a posteriori 
estimators cannot be used in the offline step. 

The global error consists of two, independent contributions: the first one is 
related to the RB approximation, the second one is related to the truncation of 
the KL expansion. 

In the greedy algorithm, we use a standard a posteriori estimation \S^ — 
^K.N n\ — ^N.K^Xm) fhe error between the finite element approximation 
^K,M ^{uK,^r{-;Yj^)) and the RB approximation S'^j^j^ := J^{uK,M,Ni-;Y^)) 
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of S'^ at a fixed truncation order K, for any realization F,^ G A^. Th is is classi- 
cal |Boy08[lNVP05llRHP] and similar to our example of Section^ see jBBM+09] 
for details. Note however that the coercivity constant of the bilinear form for 
the variational formulation 

Find u(-;2/^) e H^{V) s.t. 

[ aVu{--v'')-Vv+ [ b{-,y'')u{--y'')v^ j gv,yveH\V) (31) 



of problem (29) depends on K. To avoid the additional computation of the 
coercivity constant for each K, we impose b{x,y^) > bG{x)/2, for all x G Tb, 
and thus get a uniform lower bound for the coercivity constant. In practice, this 
imposes a limit < T < T^ax on the intensity factor in the ranges of the random 
variables Yfe, thus on the random fluctuations of the stochastic coefficient, where 
TTmax is fixed for all K £ {0, ... , K.} 

Let us now discuss the online a posteriori error estimation. As for the 
truncation error, an a posteriori estimation 15*]^ — S'^j^j-\ = \J-{u^{ ■ — 
^{uK,^r{ ■ \ Y^))\ < A%^i^{Y^) is derived in |BBM+09| . The error estima- 
tors Aj^j^(Y^) and A^^(F^), respectively for the RB approximation and 
the truncation, are eventually combined for to = 1, . . . , M to yield global error 
bounds in the Monte-Carlo estimations of the statistical outputs: \Em[{S]} ^)]- 
EM[iS]:^)]\ < AeHS^^^^n)) and |FM[(5£^,Ar)] - T/m[(5X?)]| < Ay((51?_^_^)). 
The control of the truncation error may be used to improve the performance of 
the reduced basis method. In particular, if the truncation error happens to be 
too small compared to the RB approximation error, the truncation rank /C may 
be reduced. 



3.4 Numerical results 

Our numerical simulations pre s ente d in [BBM^Pg] are performed on the steady 
heat conduction problem (23l-(24) inside the T-shaped heat sink I? C Pi U 
I?2 pictured in Fig. [IJ The heat sink comprises a 2 x 1 rectangular substrate 
(spreader) V2 = (-1, 1) x (0, 1) and a 0.5 x 4 thermal fin Vi = (-0.25,0.25) x 
(1,5) on top. The diffusion coefficient is piecewise constant, a = l-Pi + ctq 1x)2, 
where lu^ of course denotes the characteristic function of domain Vi {i — 1,2). 
The finite element approximation is computed using quadratic finite elements on 
a regular mesh, with Af = 6 882 degrees of freedom. The thermal coefficient is 
(To = 2.0. To construct the random input field b{-,Lo), we consider the covariance 
function Co\rar{b{x,uj)b{y,uj)) — (6T)^ exp(— (x — y)^/5^) for b — 0.5, T = 
0.058, and a correlation length S = 0.5. We perform its KL expansion and 
keep only the largest /C = 25 terms. We then fix G{x) = 1 and the variables 
y/c(a;), 1 < fc < /C, as independent, uniformly distributed random variables. 



This defines b{-,uj) as the right-hand side of (27). 

After computing the reduced basis offline with our RB greedy algorithm 
on a trial sample of size |Atriai| — 10 000, the global approximation error in the 
output Monte-Carlo sums EM[{S^.j\r.N)] and VM[{Sx.j\r.N)] decays very fast (in 
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Figure 2: Global error bounds for the RB approximation error and the KL 
tmncatioii error of the output expectation (top: Ae{{S^ j^))) and of the 
output variance (bottom: Ay ((S*]^^ at)))) as functions of the size = 2, . . . , 14 
of the reduced basis, at different truncation orders K = 5, 10, 15, 20. 



22 



fact, exponentially) with the size TV = 1, . . . , 14 of the reduced basis, see Fig. [2] 
with K = 20. Note that M = 10 000 for the Monte-Carlo sums. We would also 
like to mention that these reduced bases have actually been obtained letting 
varying not only the parameter , but also additional parameters (namely 
the diffusion coefficient a and the mean b of the Biot number) but this does 
not influence qualitatively the results presented here, and we omit this technical 
issue for simplicity (see |BBM+09) for more details). 

It is observed that the global approximation error for truncated problems 
at a fixed order K and for various N (the size of the reduced basis) is quickly 
dominated by the truncation error. More precisely, beyond a critical value 
N > Ncrit{K), where iVcrit(-fl^) is increasing with K, the global approximation 
error becomes constant. Notice that the approximation error is estimated as 
usual by a posteriori estimation techniques. 

When /C is infinite (or finite but huge) , the control of the KL truncation error 
may be difficult. This is a general issue for problems involving a decomposition 
of the stochastic coefhcient. Our RB approach is still efficient in some regimes 
with large K, but not all. In particular, a fast decay of the ranges of the 
parameters {yk)i<k<K facilitates the exploration of by the greedy algorithm, 
which allows in return to treat large K when the eigenvalues Xk decay sufficiently 
fast with k. 

In [BBM+09| . we have decreased the correlation length to 6 = 0.2 and could 
treat up to A" = 45 parameters, obtaining the results in a total computational 
time still fifty times as short as for direct finite element computations. 



4 Variance Reduction using an RB approach 

In this section, we present a variance reduction technique based upon an RB 
approach, which has been proposed recently in |BL09) . In short, the RB ap- 
proximation is used as a control variate to reduce the variance of the original 
Monte-Carlo calculations. 



4.1 Setting of the problem 

Suppose we need to compute repeatedly, for many values of the parameter A G A, 
the Monte-Carlo approximation (using an empirical mean) of the expectation 
E(Z'^) of a functional 

Z^=g\X^)- r f\s,X^)ds (32) 

"'0 

of the solutions [X^,t e [0,r]) to the Stochastic Differential Equation (SDE) 

X^=x+f b^{s,X^)ds+ f a^{s,X^)dB„ (33) 
Jo Jo 

where (Bt eR'^,te [0,T]) is a d-dimcnsional standard Brownian motion. The 
parameter A parametrizes the functions g^, /^, fe^ and cr^. In (33), we assume 
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6^ and allow for the Ito processes {X^ €W^,t€ [0, T]) to be well defined, for 
every A G A. Notice that we have supplied the equation with the deterministic 
initial condition = x ^W^. In addition, and are also assumed smooth, 



such that Z G L (il). Recall that a symbolic concise notation for (33) is 



dX^ = b^{t, X^) dt + a^{t, X^) dBt with X^ = X. 

Such parametrized problems are encountered in numerous applications, such 
as the calibration of the volatility in finance, or the molecular simulation of 
Brownian particles in materials science. For the applications in finance, E(Z'*') 
is typically the price of an European option in the Black-Scholes model, and 
A enters the diffusion term (the latter being called the volatility in this con- 
text). The calibration of the volatility consists in optimizing A so that the 
prices observed on the market are close to the prices predicted by the model. 
Any optimization procedure requires the evaluation of E(Z^) for many values 
of A. On the other hand, the typical application we have in mind in materi- 
als science is related to polymeric fluids modelling. There, 'Ei{Z^) is a stress 
tensor which enters the classical momentum conservation equation on velocity 
and pressure, and X^ is a vector describing the configuration of the polymer 
chain, which evolves according to an overdamped Langevin equation, namely a 



stochastic differential equation such as (33). In this context, A is typically the 



gradient of the velocity field surrounding the polymer chain at a given point in 
the fluid domain. The parameter A enters the drift coefficient b^. The com- 
putation of the stress tensor has to be performed for each time step, and for 
many points in the fluid domain, which again defined a many-query context, 
well adapted to the RB approach. For more details on these two applications, 
we refer to jBL09| . 



We consider the general form (|32|)-(33) of the problem and as output the 



Monte-Carlo estimation E]vi[(-^m)] = i/X^m^i-^m parametrized by A e A, 
where we recall {Z^) denotes i.i.d. random variables with the same law as 
Z^ . These random variables are build in practice by considering a collection 



of realizations of (33), each one driven by a Brownian motion independent 
from the others. In view of the Central Limit Theorem, the rate at which 
the Monte-Carlo approximation Em [(.^4)] approaches its hmit E(Z'^) is given 
by the prefactor being proportional to the variance oi Z^ . A standard 

approach for reducing the amount of computations is therefore variance re- 
duction |Aron41 IM()951 IOvdBH97l IBPMI IHDH641 IMTOBj . We focus on one 
particular variance reduction technique: the control variate method. It consists 
in introducing a so called control variate G L^(r2), assumed centered here 
for simplicity: 

E(y^) = 0, 

and in considering the equality: 

E(Z^) = E{Z^ -Y^). 

The expectation E(Z'^ — y'^) is approximated by Monte-Carlo estimations Em[(.Z',^ 
Yj^^)] which hopefully have, for a well chosen Y^, a smaller statistical error than 
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direct Monte-Carlo estimations EM[(^m)] of 'E{Z^). More precisely, is ex- 
pected to be chosen so that Var(Z^) > Var(Z^ - Y^). The central hmit 
theorem yields 

M 

Em[(^,^, - Y„\)] := ^ E (^,^, - Y^) EiZ^ Y^), (34) 

ni—1 

where the error is controlled by confidence intervals, in turns functions of the 
variance of the random variable at hand. The empirical variance 

VarM ((^;^ - r,^)) := i^n - Y,^ - Em{{Z^ - Y^))f (35) 

which, as M — ^ oo, converges to Var(Z^), yields a computable error bound. 
The Central Limit Theorem indeed states that: for all a > 0, 

P (^|E(Z^ - Y^) Em {[Zl - r,^J) I < ,y VarM {{Z^ - y^) ) ^ e_ 

(36) 

Evaluating the empirical variance (35 1 is therefore an ingredient in Monte-Carlo 
computations, similar to what a 'posteriori estimates are for a deterministic 
problem. 

Of course, the ideal control variate is, VA e A: 

Y^ = Z^ - E(Z^) , (37) 

since then, Var {Z^ — F^) = 0. This is however not a practical control vari- 
ate since E(Z'^) itself, the quantity we are trying to evaluate, is necessary to 
compute (37 1. Ito calculus shows that the optimal control variate (37) also 
writes: 



Y^ = [ Vu^is,X^) ■ a^{s,X^)dB, 
Jq 



(38) 



where u^{t,y) G {[Q,T],C'^{R'^)) satisfies the backward Kolmo gorov equa- 
tion: 

1 u\T,-)=g^{-). 



Even using this reformulation, the choice (37 1 is impractical since solving the 
partial differential equation (391 is at least as difficult as computing E(Z^). We 
will however explain now that both "impractical" approaches above may give 
birth to a practical variance reduction method, when they are combined with a 
RB type approximation. 

Loosely speaking, the idea consists in: (i) in the offiine stage, compute fine 
approximations of E(Z'^) or respectively for some appropriate values of A, in 
order to obtain fine approximations of the optimal control variate Y^ (at those 
values) and (ii) in the online stage, for a new parameter A, use as a control 
variate the best linear combination of the variables built offline. 
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4.2 Two algorithms for variance reduction by the RB ap- 
proach 

Using suitable time discretization methods [KPOO' , realizations of the stochastic 



process ( 33 ) and the corresponding functional ( 32 ) can be computed for any 
A € A, as precisely as needed. Leaving aside all technicalities related to time 
discretization, we thus focus on the Monte Carlo discretization. 
We construct two algorithms, which can be outlined as follows. 



Algorithm 1 (based on formulation (37)) 



Offline stage: Build an appropriate set of values {Ai, . . . , Ajv} and, con- 
currently, for each A G {Ai, . . . , Xn} compute an accurate approximation 
-f^Miarge, [(^m)] ^f E(Z'^) (for & vcry large number Miargo of realizations). 
At the end of the offline step, accurate approximations 

of the optimal control variate are at hand. The set of values {Ai, . . . , Aat} 
is chosen in order to ensure the maximal variance reduction in the forth- 
coming online computations (see below for more details). 

Online stage: For any A G A, compute a control variate Y^ for the Monte- 
Carlo estimation of E(Z^) as a linear combination of 

(y, = f^-)i<.<w. 



Algorithm 2 (based on formulation (38)) 



Offline stage: Build an appropriate set of values {Ai, . . . , Aat} and, con- 
currently, for each A G {Ai, . . . , Xn}, compute an accurate approximation 
tt'^ of u^, by solving the partial differential equation (|39|. The set of values 
{Ai, . . . , A^r} is chosen in order to ensure the maximal variance reduction 
in the forthcoming online computations (see below for more details). 

Online stage: For any A G A, compute a control variate Y^ for the Monte- 
Carlo estimation of E(Z^) as a linear combination of 



l<i<N 



In both algorithms, we denote by Y^ the control variate built online as a linear 
combinations of the Yi^s, the lowerscript index TV emphasizing that the ap- 
proximation is computed on a basis with N elements. An important practical 
ingredient in both algorithms is to use for the computation of Z^ the exact same 
Brownian motions as those used to build the control variates. 

The construction of set of values A G {Ai, . . . , Ajv} in the offline stage of 
both algorithms is done using a greedy algorithm similar to those considered in 
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the preceeding sections. The only diflFerence is that the error estimator used is 
the empirical variance. Before entering that, we need to make precise how the 
linear combinations are built online, since this linear combination construction 
is also used offline to choose the A^'s. 

The online stage of both algorithms 1 and 2 follow the same line: for a given 
parameter value A € A, a control variate for Z\ is built as an appropriate 
linear combination of the control variates (yi)i<i<jv (obtained from the offline 
computations). The criterium used to select this appropriate combination is 
based on a minimization of the variance: 

N 

Y^ = Va:F„, (40) 



where 



(a*)i<„<Ar = arg min Var \ Z -} ^a^Yn ) . (41) 

{an)l<n<N^ 



N 
n=l 



In practice the variance in (41 1 is of course replaced by its empirical approx- 
imation VarM,„^u- Notice that we have an error estimate of the Monte Carlo 
approximation by considering VarAf^„^,i (^Z^ — Y^^ . It is easy to check that 

the least squares problem (41 1 is computationally inexpensive to solve since it 
amounts to solving a linear N x N system, with N small. More precisely, this 
linear system writes: 

where a* here denotes the vector with components a* , Cm^^^u is a matrix with 
(i, j)-th entry 

CovAf3„,,i (Yi 
and fcAfsmaii is a vecteur with j-th component 

C0VA4„,^„(Z4,yj-,m) 

where for two collections of random variables Um and Vm, 

1 

C0VA/(J7m,Ki) = E ^rnVm - 

In summary, the computational complexity of one online evaluation is the sum of 
the computational cost of the construction of &Ms„aii (wich scales like iVMsmaii), 
and of the resolution of the linear system (which scales like N"^ for Algorithm 1 
since the SVD decomposion of Cm^^^u be precomputed offline, and scales 
like for Algorithm 2, since the whole matrix Cm^^^h has to be recomputed 
for each new value of A) . 

The greedy algorithms used in the offline stages follow the same line as in 
the classical RB approach. More precisely, for Algorithm 1, the offline stage 
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writes: Let Ai S Atrial be already chosen and compute EM,„g„(^'^0- Then, for 
i = 1, . . . , N — 1, for all A € Atrial, compute Y/" and inexpensive approximations: 

E,iX) := Em.„„„(^^ - f,^) for E(Z^) , 
e,(A) := VarM.„.„ (^^ - f,^) for Var (z^ - K,^) . 
Select Ai+i e argmax {ei(A)}, and compute EM,„gc (^'^'"''O- 

AeAt,iai\{Aj J = l,...,i} 

In practice, the number N is determined such that eAr(AAr+i) < e, for a given 
threshold e. The greedy procedure for Algorithm 2 is similar. 

4.3 Reduced-Basis ingredients 

The algorithms presented above to build a control variate using a reduced basis 
share many features with the classical RB approach. The approach follows a 
two-stage offline / online strategy. The reduced basis is built using snapshots 
(namely solutions for well chosen values of the parameters). An inexpensive 
error estimator is used both in the offline stage to build the reduced basis in the 
greedy algorithm, and in the online stage to check that the variance reduction 
is correct for new values of the parameters. The construction of the linear 
combinations for the control variates is based on a minimization principle, which 
is reminiscent of the Galerkin procedure ^ . 

The practical efficiency observed on specific examples is similar for the two 
algorithms. They both satisfactorily reduce variance. Compared to the plain 
Monte Carlo method without variance reduction, the variance is divided at 
least by a factor 10^, and typically by a factor 10^. Algorithm 2 appears to 
be computationally much more demanding than Algorithm 1 and less general, 
since it requires the computation (and the storage) of an approximation of the 
solution to the backward Kolmogorov equation p9| ) for a few values of the 
parameter. In particular. Algorithm 2 seems impractical for high dimensional 
problems {X^ £ M'' with d large). On the other hand. Algorithm 2 seems to 
be more robust with respect to the choice of Atrial : it yields good variance 
reduction even for large variations of the parameter A, in the online stage. We 
refer to jBLOQj for more details. 

Notice also that Algorithm 1 is not restricted to a random variable that is 
defined as a functional of a solution to a SDK. The approach can be generalized 
to any parametrized random variables, as long as there is a natural method 
to generate correlated samples for various values of the parameter. A natural 
setting for such a situation is the computation of a quantity K{g'^{X)) for a 
random variable X with given arbitrary law, independent of the parameter A. 
In such a situation, it is easy to generate correlated samples by using the same 
realizations of the random variable X for various values of the parameter A. 
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4.4 Numerical Results 

The numerical results shown on Figure [s] are taken from |BL09| and relate to the 
second application mentioned in the introduction, namely multiscale models for 
polymeric fluids (see |LL07j for a general introduction). In this context, the non- 
Newtonian stress tensor is defined by the Kramers formula as an expectation 
E{Z^) of the random variable: 

=X^(g) F(X^) , (42) 

where is a vector modelling the conformation of the polymer chain. The 
latter evolves according to an overdamped Langevin equation: 

dX^ = (A X^ - F{X^)) dt + dBf (43) 



Equation (43) holds at each position of the fluid domain, the parameter A G 
T^d-xd (^fi — 2 or 3) being the local instantaneous value of the velocity gradient 
field at the position considered. The evolution of the "end-to-end vector" X^ 
is governed by three forces: a hydrodynamic force XX^, Brownian collisions 
Bt against the solvent molecules, and an entropic force F{X^) specific to the 
polymer molecule. Typically, this entropic force reads either F{X^) — X^ (for 

the Hookean dumbbells), or F{Xi') = i_^x*>-\2/b C^*-*^ Finitely-Extensible 

Nonlinear Elastic (FENE) dumbells, assuming \X^\ < Vb). 

The numerical simulations of the flow evolution of a polymeric fluid using 
such a model typically consist, on many successive time slots [nT, [n + 1)T], of 



two steps: (i) the computation of (431, for a given gradient velocity field A, at 
many points of the fluid domain (think of the nodes of a flnite element mesh) 
and (ii) the computation of a new velocity gradient field in the fluid domain, 
for a given value of the non-Newtonian stress tensor, by solving the classical 
momentum and mass conservation equations, we omit here for brevity. Thus, 
E(Z'^) has to be computed for many values A corresponding to many spatial 
positions and many possible velocity fields at each such positions in the fiuid 
domain. 

In the numerical simulations of Figure [3j the SDE (43) for FENE dumbbells 
when c? = 2 is discretized with the Euler-Maruyama scheme using 100 iterations 
with a constant time step At = 10~^ starting from a deterministic initial con- 
dition X = (1, 1). Refiecting boundary conditions are imposed on the boundary 
of the ball with radius Vb. For & = 16 and | Atrial | = 100 trial parameter values 
randomly chosen in the cubic range A = [— 1, 1]^ (the traceless matrix A has en- 
tries (All = ^-^22, Ai2, A21)), a greedy algorithm is used to incrementally select 



= 20 parameter values after solving | Atrial | = 100 least-squares problems (41 ) 
(with Msinaii = 1000) at each step of the greedy algorithm (one for each of the 
trial parameter values A £ Atrial)- Then, the A'' = 20 selected parameter values 
are used online for variance reduction of a test sample of |Atcst| = 1000 random 
parameter values. 

The variance reduction obtained online by Algorithm 1 with Afiargc — lOOMsmaii 
is very interesting, of about 4 orders of magnitude. For the Algorithm 2, we 
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use the exact solution to the Kolmogorov backward equation for Hookean 
dumbells as an approximation to u"^ solution to (39). This also yields satisfy- 
ing variance reduction though apparently not as good as in Algorithm 1. As 
mentioned above, Algorithm 2 is computationally more demanding but seems 
to be slightly more robust than Algorithm 1 (namely when some online sample 
test Atestwidc Uniformly distributed in [—2, 2]'^ extrapolates the trial sample used 
offline, see Fig. |3|. 
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Figure 3: Algorithm 1 (left) and 2 (right) for FENE model with b = 16. The 
X-axis is the size N of the reduced basis. We represent the minimum +, mean 
X and maximum o of VarM[^'*' — Y^]/'Em[Z^ — Y^]'^ over online test samples 
Atost C A (top) and Atestwide 3 A (bottom) of parameters. 

Our numerical tests, although preliminary, already show that the reiterated 
computations of parametrized Monte-Carlo estimations seem to be a promising 
opportunity of applications for RB approaches. More generally, even if RB 
approaches may not be accurate enough for some applications, they may be 
seen as good methods to obtain first estimates, which can then be used to 
construct more refined approximations (using variance reduction as mentioned 
here, or maybe preconditionning based on the coarse-grained RB model). This 
is perhaps the most important conclusion of the work described in this section. 
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5 Perspectives 



The standard RB method has proved numerically efhcient and reliable at reduc- 
ing the cost of computations for the approximation of solutions to parametrized 
boundary value problems in numerous benchmark many-query frameworks. These 
accomplishments claim for a wider use of the RB method in more realistic set- 
tings, and even suggest that some RB ideas could still be extended in numerous 
many-query frameworks yet largely unexplored, including the stochastic context. 
The success of the RB approach in parametrized boundary value problems is 
only understood precisely from a mathematical perspective in a few very simple 
cases. This should motivate further theoretical investigations. 

In this section, we discuss various tracks for the development of reduced 
basis techniques, both from a methodological viewpoint and in terms of possible 
applications, with a focus on the stochatic context presented above. 

5.1 A posteriori estimation in the stochastic context 

We already emphasized that a crucial ingredient in the RB approach is an ac- 
curate and fast a posteriori estimator for the approximation error between two 
levels of discretization (the initially discretized, non-reduced one and the re- 
duced one). Therefore, before everything, the future developments of the RB 
method should definitely concentrate on improving the a posteriori estimators. 
In particular, for the new contexts of application that are stochastic, there seems 
to remain some room for a yet better understanding of the a posteriori error 
estimation. More precisely, the best way to evaluate the reduction error when 
the Galerkin approximations (used by deterministic applications) are replaced 
with Monte-Carlo approximations is still unclear. For a first application of 
the RB ideas to stochastic applications, we have used confidence intervals as 
a probabilistic measure of the Monte-Carlo approximation error. These confi- 
dence intervals are only reliable in the limit of infinitely many realizations of the 
random variables. But there are other possibilities, like using non- asymptotic 
upper-bounds for the error which hold whatever the number of realizations (us- 
ing for example Chebyshev inequalities or Berry- Esseen type bounds). In ad- 
dition, the numerical evaluation of the variance is not obvious either. Until 
now, we have used Monte-Carlo estimators, but there exist other possibilities 
too which could be faster or more accurate and should thus be tested. Finally, 
another idea related to the method presented in Section [4] would be to mimick 
the usual RB approach, by considering that the reference result is the one ob- 
tained with Miarge realizations, and to develop a posteriori error bounds with 
respect to this reference solution (using for example conditional expectations 
with respect to the Mjarge realizations). 

5.2 AfRne decompositions and the stochastic context 

As explained above, the RB approach is to date only efficient at yielding compu- 
tational reductions in the context of affine parametrization. However, as shown 
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in the previous sections, a many-query parametrized framework is not necessar- 
ily parametrized in an affine way. So one may have to pretreat the problem in 
order to transform it as the limit of a sequence of affinely-parametrized problems. 
It would thus be interesting to derive rapidly convergent affine approximations 
for non-affine problems. For instance, the Karhunen-Loeve decomposition used 
to pretreat a random field entering a partial differential equation as a coefficient 
may converge too slowly for an efficient use of the RB method applied to trun- 
cated decompositions, in the context of random fields with small correlation 
lengths. One should then look for other possible affine representation of the 
random variations in the input coefficient. Now, there are many possible tracks 
to solve this problem, like projecting the random field on a well-chosen basis 
for the oscillation modes of the coefficient for instance (many possible bases 
may exist for the realizations of the random field, depending on its regularity), 
interpolating (recall the so-called empirical interpolation with magic points), 
etc. 

5.3 Application to Bayesian statistics 

A context where the RB ideas could be applied is Bayesian statistics where the 
many-query parametrized framework is naturally encountered. We now make 
this more precise by presenting a specific example that would be well-suited for 
the application of Algorithm 1 in Section |4] 

Let us consider, for given values of a parameter fii, an ensemble of observa- 
tions (2;f^)]^<j<^f'i . Following the Bayesian framework, a stochastic model is 
proposed to model the observation: the quantities (2;f^)i<i<Ar'^i are supposed 
to form a set of independent and identically distributed samples following a given 
distribution parametrized by another set of parameters 112 (think for example 
of a mixture of Gaussians, ^2 being then the triplets of weights, means and 
variances of each Gaussians) . The Bayesian approach then consists in postulat- 
ing a so-called prior distribution (with a probability density function denoted 
Prior(/i2) below) on the parameters ^2, and to compute the so-called posterior 
distribution, namely the distribution of fj,2 given the observations (with a prob- 
ability density function denoted Il[^2\{x^^)l<:i<:J~^^'l ) below). Of course, the 
posterior distribution for /i2 is expected to depend on fii: for each fii, the aim 
is thus to sample the probability measure n(/i2|(2;f^)i<i<7vfi ) d/i2i with 

n(/^2|«^)i<.<A,^-i^J = (Z^^)-in(«^)i<,<^.^.^jM2) Prior(/i2) 

where Z'^^ is the normalization constant, and n((a;f^)i<i<Arf'i |/^2) is the so- 
called likelihood function, namely the probability density function of the obser- 
vations given the datas. One possible technique to sample the posterior distri- 
bution consists in drawing samples according to the prior distribution, and to 
weight each of them using the likelihood function, which depends on fii. With 
such a sampling technique, it is easy to draw correlated samples for various 
values of /ii . Following Algorithm 1 in Section |4j it would thus be possible to 
build a reduced basis based on the sampling of the posterior distribution for 
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some selected values of (offline stage), in order to reduce the variance for 
the sampling of the posterior distribution for other generic values of /ii (online 
stage). 

5.4 Relation to functional quantization 

One computationally demanding stochastic context that defines a many-query 
framework is the approximation of the solution to a parametrized stochastic 
differential equation, for many values of the parameter. We already mentioned 
applications in finance and rheology in Section 4, where a variance reduction 
technique based on RB was proposed. Another idea consists in first computing 
precise discretizations of a few processes at some well-chosen parameter values, 
and then to use them for a faster computation of an approximation of the 
processes for other values of the parameter. 

Functional quantization is an approach that has independently been devel- 
oped along this line, see for instance |LuPa021 IPPOSj . The idea of quantization 
is to approximate a square-integrable random variable with values in a Hilbert 
space by a random variable that takes a finite number of values, in an optimal 
way. In its simplest form, quantization deals with Gaussian random variables 
with values in E'^, but it can also be applied to Gaussian processes. The nu- 
merical approach developped in |PS08| to solve stochastic differential equations 
is to first quantize the Brownian motion, and then to solve a collection of ordi- 
nary differential equations in order to recover approximations of the solutions 
to the stochastic differential equations as linear combinations of the ordinary 
differential equations solutions. Clearly, this approach for the discretization of 
stochastic differential equations has intimate connection with a RB approach. 
In particular, the computations are split into two parts: an offline step, which 
is computationally expensive, to quantize the Brownian motion, and then an 
online step to solve ordinary differential equations rather than stochastic differ- 
ential equations. 

In a setting where the stochastic differential equations are parametrized, a 
natural similar idea would be to quantize the solution to the stochastic differen- 
tial equations for a few values of the parameter, and next to build the solution 
to the stochastic differential equation for another value of the parameter as a 
linear combination of these precomputed solutions. 
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